######################################################################
#Replication file for "The Impact of China's AIIB on the World Bank"
#Authors: Jing Qian, James Vreeland, Jianzhi Zhao
#Codes to replicate Figure A.10
######################################################################
rm(list = ls())

#To install the bpCausal package (ver 0.0.1)
# install.packages('devtools', repos = 'http://cran.us.r-project.org') # if not already installed
# devtools::install_github('liulch/bpCausal')

#Load bpCausal package
library(bpCausal)

#Other packages used
library(tidyverse)
library(RColorBrewer)

#Setwd
setwd("C:/Users/qianj/Dropbox (Princeton)/AIIB_WB_Replication")

#Load custom functions
#Note: Given the number of tests performed in the appendix, 
#several custom functions are created to implement estimation methods in a more efficient way.
#These functions merely create a wrap-up of existing functions, like bpCausal and gsynth,
#in order to avoid typing parameters repeatedly.

source("code/custom_functions.R")

###################################
#Estimation for Figure A.10
###################################
#----------------
#(1) Load data
#----------------
load("data/data_main.RData")

#----------------
#(2) Setup
#----------------
df.use = df.main
D.use = "aiib_founder_2016"
covs.all = c("gdppc_log_lag",
             "population_log_lag",
             "election_either_lag",
             "fdi_gdp_lag",
             "debt_gni_lag",
             "oda_gni_lag",
             "polity2_lag",
             "unsc_lag",
             "IdealPointDistance_lag")
X.use = Z.use = A.use = covs.all
threshold = seq(50, 100, 10)

#Container
result.a10 = list()

#----------------
#(3) Estimation
#~ 35 minutes
#----------------
for (i in 1:length(threshold)){
  dv.list = paste0(c("hard", "soft"),
                   "_project_count_all_",
                   threshold[i])
  result.a10[[i]] = bpcausal.group(dv.list = dv.list,
                                   df.use = df.use,
                                   D.use = D.use,
                                   X.use = X.use,
                                   Z.use = Z.use,
                                   A.use = A.use)
  print(i)
}


###################################
#Produce Figure A.10
###################################
#-------------
#(1) Prepare results
#-------------
#Summarize
result.threshold.use = lapply(result.a10,
                              function(x) lapply(x, effSummary))

#-------------
#(2) Parameters
#-------------
#Parameters
n = 2
xlim = c(-2, 2)
cex = 2
xtext = 0.65
colors = brewer.pal(n = 6,
                    name = "Dark2")
k = 6
gap = 1/(k+2)

#------------------
#(3) Plot
#------------------
#----------
#Setup
png(filename = "figure/Figure_A10.png",
    height = 800,
    width = 1200)

par(mar = c(5, 15, 1, 1))

#----------
#Empty plot
plot(1,
     type = "n",
     xlab = "Estimated Average Treatment Effect on the Treated",
     ylab = "",
     yaxt = "n",
     xaxt = "n",
     xlim = xlim,
     ylim = c(n + 0.5, 0.5),
     cex.lab = 2)
box()

#------------------
#Names
#------------------
percentage = seq(50, 100, 10)
names.hard = paste0("hard_project_count_all_",
                    percentage)
names.soft = names.hard %>%
  gsub("hard", "soft", .)

#------------------
#Hard
#------------------
for (i in 1:length(names.hard)){
  points(x = result.threshold.use[[i]][[names.hard[i]]]$est.avg[1],
         y = 1 + (i-3)*gap,
         col = colors[i],
         pch = 19,
         cex = cex)
  arrows(x0 = result.threshold.use[[i]][[names.hard[i]]]$est.avg[2],
         x1 = result.threshold.use[[i]][[names.hard[i]]]$est.avg[3],
         y0 = 1 + (i-3)*gap,
         y1 = 1 + (i-3)*gap,
         length = 0.05,
         angle = 90,
         code = 3,
         col = colors[i],
         cex = cex,
         lwd = cex)
}


#------------------
#Soft
#------------------
for (i in 1:length(names.soft)){
  points(x = result.threshold.use[[i]][[names.soft[i]]]$est.avg[1],
         y = 2 + (i-3)*gap,
         col = colors[i],
         pch = 19,
         cex = cex)
  arrows(x0 = result.threshold.use[[i]][[names.soft[i]]]$est.avg[2],
         x1 = result.threshold.use[[i]][[names.soft[i]]]$est.avg[3],
         y0 = 2 + (i-3)*gap,
         y1 = 2 + (i-3)*gap,
         length = 0.05,
         angle = 90,
         code = 3,
         col = colors[i],
         cex = cex,
         lwd = cex)
}
#------------------
#Other elements
#------------------
#Abline
abline(v = 0,
       lty = "dashed",
       lwd = cex,
       col = "gray")
#Legend
legend("topright",
       legend = c(paste0(">=", 
                         percentage[-6],
                         "%"),
                  "=100%"),
       cex = cex,
       col = colors,
       lty = c(1, 1),
       bty = "n",
       lwd = cex)

#Axis
axis(side = 1,
     at = seq(min(xlim),
              max(xlim),
              1),
     cex.axis = cex)

text(x = min(xlim) - xtext,
     y = (1:2) + 0.1,
     xpd = NA,
     cex = cex,
     labels = c("Infrastructure",
                "Non-Infrastructure"))
#----------
#Close
dev.off()